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Abstract 

A one-dimensional branching-coalescing model is considered on a chain of length 
L with reflecting boundaries. We study the phase transitions of this model in a 
canonical ensemble by using the Yang-Lee description of the non-equilibrium phase 
transitions. Numerical study of the canonical partition function zeros reveals two 
second-order phase transitions in the system. Both transition points are determined 
by the density of the particles on the chain. In some regions the density profile of 
the particles has a shock structure. 

Key words: Yang-Lee Theory, Matrix Product Formalism, Shock 
PACS: 05.40.-a, 05.70.Fh, 02.50.Ey 



One-dimensional driven lattice gases are models of particles which diffuse, 
merge and separate with certain probabilities on a lattice with open, periodic 
or reflecting boundaries. In the open boundaries case the particles are allowed 
to enter or leave the system from both ends or only one end of the chain. In 
the reflecting boundaries or periodic boundary cases the number of particles 
will be a conserved quantity provided that no other reactions other than the 
diffusion of particles take place. In the stationary state, these models exhibit 
a variety of interesting properties such as non-equilibrium phase transitions 
and spontaneous symmetry breaking which cannot be found in equilibrium 
models (see [1] and references therein). Different applications are also found 
for such models which include the kinetics of biopolymerization [2] and traffic 
flow modelling [3]. These models have also allowed the study of shocks i.e. 
discontinuities in the density of particles over a microscopic region. Over the 
last decade people have studied the shocks in one-dimensional driven-diffusive 
models with open and periodic boundary conditions. A prominent example of 
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such models with open boundary is the Asymmetric Simple Exclusion Process 
(ASEP) in which the particles enter the system from the left boundary, diffuse 
in the bulk and leave the chain from the right boundary with certain rates 
[4]. For specific tuning of injection and extraction rates a shock might appear 
in the system which moves with a constant velocity towards the boundaries. 
The ASEP has also been studied on a ring in the presence of a second class 
particle called the impurity [5,6,7]. In this case the impurity will track the 
shock front with a constant velocity which is determined by the reaction rates 
of the model. The shocks in the models with reflecting boundaries have not 
been studied yet. 

In the present letter we study the phase transitions in a one-dimensional 
branching-coalescing model with reflecting boundaries in which the particles 
diffuse, coagulate and decoagulate on a lattice of length L. The reaction rules 
are specifically as follows: 

Diffusion to the left: + A -> A + with rate q 

Diffusion to the right: A + -> + A with rate q' 1 

Coalescence to the left: A + A — > A + with rate q 

Coalescence to the right: A + A — > + A with rate q' 1 

Branching to the left: + A -> A + A with rate Aq 

Branching to the right: A + — > A + A with rate Aq' 1 



(1) 



in which A and stand for the presence of a particle and a hole respectively. 
It is assumed that there is no injection or extraction of particles from the 
boundaries. We will also assume that the number of particles on the chain 
is a conserved quantity. This model was first introduced and treated in con- 
tinuum approximation in [8,9,10,11]. It was then studied using the Empty 
Interval Method (EIM) in [12] . In this formalism the physical quantities such 
as the density of particles are calculated from the probabilities to find empty 
intervals of arbitrary length. Later this model was studied using so-called the 
Matrix Product Formalism (MPF) [13]. According to this formalism the sta- 
tionary probability distribution function of the system is written in terms of 
the products of non-commuting operators E and D and the vectors \V) and 
(W\ as follows 

P(n, -~,t l ) = ±{w\ flinD + (i - r t )E)\v). (2) 

Z L i=1 



Each site of the lattice is occupied by a particle (t{ = 1) or is empty (r, = 0). 
The factor Zl in (2) is a normalization factor. The operators D and E stand 
for the presence of particles and holes respectively and besides the vectors \ V) 
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and (W\ should satisfy the following quadratic algebra [13] 



[E, E] = 

ED-ED = q(l + A) ED - \DE - \D 2 
DE — DE = —qED + ^DE - qD 2 
DD- DD = -qAED - ^DE + (q + ±)L> 2 
(W\E= (W\D = 
E\V) = D\V) = 0. 



The operators D and E are auxiliary operators and do not enter into calculat- 
ing (2). Having a representation for the quadratic algebra (3) one can easily 
compute the steady state weights of any configuration of the system using 
(2). It has been shown that (3) has a four-dimensional representation [13]. For 
q 2 7^ 1 + A we have 
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in which a and b are arbitrary constants. The matrix representations for D 
and E are also given in [13]. Both EIM and MPA approaches showed that the 
model has two different phases: a low-density phase for q 2 > 1 + A and a high- 
density phase for q 2 < 1 + A and a phase transition takes place at the critical 
point where q 2 = 1 + A. In the low-density and the hight-density regions the 
density profile of the particles on the chain is an exponential function while 
on the coexistence line it changes linearly along the lattice. 
Recently this model has been studied under the open boundary conditions on 
a specific manifold of the parameters of the system [20] . It is shown that if the 
particles are injected and extracted from the left boundary with the rates a 
and f3 respectively then the model has the same phase structure provided that 
a = (q' 1 — q + P)A. In the latter case the operators D and E and also the vec- 
tors (W\ and \ V) have two-dimensional representations. The only difference is 
that for the reflecting boundary conditions the system involves three different 
length scales while for the open boundary conditions it is characterized by 
one length scale. Moreover, it has been shown that in the open boundary case 
the probability distribution function of the system can be written in terms of 
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superposition of Bernoulli shock measures [19]. For the open boundary case if 
we fix the density of particles, for example by working in a canonical ensem- 
ble, we can see the real shock structures in the density profile of particles. In 
this case the system will also have two different phases: a low-density phase 
and a jammed-phase where the shocks evolve in the system. These phases are 
specified by the density of particles and are separated by a second-order phase 
transition [21]. 

A natural question that might arise is whether or not we can see the shocks in 
our branching-coalescing model defined by (1) with the reflecting boundaries. 
To answer this question we will study the model with reflecting boundaries 
in a canonical ensemble where the number of particles on the chain is equal 
to M so that the density of particles p — ^ remains constant. We will then 
investigate the phase transitions and the density profile of particles on the 
chain. Recently it has been shown that the classical Yang-Lee theory [14,15] 
can be applied to the one-dimensional out-of-equilibrium systems in order to 
study the possible phase transitions of these models [16,17,18,20,21]. Accord- 
ing to this theory in the thermodynamic limit, the zeros of the canonical or 
grand canonical partition function, as a function of an intensive variable of 
the system, might approach the real positive axis of that parameter at one 
or more points. Depending on how these zeros approach the real positive axis 
the system might have one or more phase transition of different orders. If the 
zeros intersect the real positive axis at a critical point at an angle then n 
will be the order of phase transition at that point [16]. 
Let us define the canonical partition of our model using the MPF as follows 

Z LM = E S(M-j2T t )(W\f[(T t D + (l-T t )E)\V) (5) 
{n=0,l} i=l i=l 

in which M and L are the number of particles and the length of the system 
respectively and 8(- ■ •) is the ordinary Kronecker delta function 5 Xj0 . Using the 
matrix representations D, E, (W\ and \ V) given by (4) we have been able to 
calculate the canonical partition of this model (5) and its zeros numerically. 
One can use MATHEMATICA to calculate (W\(E + xD) L \V) for arbitrary 
q and A and finite L in which x is a free parameter. The result will be a 
polynomial of x. The coefficient x M in this polynomial gives the canonical 
partition function of the system. Formally we can write 

Z LM = Coefficient[(W\(E + xD) L \V),M] (6) 

in which Coef ficient[Expr, n] gives the coefficient of x n in the polynomial 
Expr. In Fig. 1 we have plotted the numerical estimates for the zeros of Zl,m 
obtained from (6) on the complex-g plane for L = 80, M = 42. The canonical 
partition function (6) has 4(L — M) zeros in the complex-g plane. We have 
found that for large L and M the locations of these zeros are not sensitive 
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to the value of A. We have also calculated the numerical estimates for the 
roots of (6) as a function of A for fixed values of q. It turns out that (6), as a 
function of A, does not have any positive root; therefore, we expect that the 
phase transition points do not depend on A. As can be seen in Fig. 1 the zeros 
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Fig. 1. Plot of the numerical estimates for the canonical partition function zeros 
obtained from (6) for L = 80 and M = 42. 

lie on two different curves and accumulate towards two different points on the 
positive real-g axis. By extrapolating the real part of the nearest roots to 
the positive real-g axis for large L and M, we have found that the transition 
points are q c = ^== (1 < q c < oo) and q' c = y/\ — p (0 < q' c < 1). As 
p — > the two curves lie on each other and we will find only one transition 
point at q c — q' c — 1. It appears also that the zeros on both curves approach 
the real-g axis at an angle | (the smaller angle). This predicts two second- 
order phase transitions at q c and q' c . The reason that the system has two 
phase transitions can easily be understood. The parameter q determines the 
asymmetry of the system and for any q the system is invariant under the 
following transformations 



q >q 

i — >L-i + l. (7) 

Therefore, one can expect to distinguish two critical points which are related 
according to the symmetry of the system. 

Let us now study the density profile of particles on the chain p(i) in each 
phase. The density of particles at site i is defined as 

T,c^{ r i, ■■■,r L ) 



in which C is any configuration of the system with fixed number of particles 
M and P(t 1: • • • ,t l ) is given by (2). It can be verified that the density profile 
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of the particles p(i) can be written as 

(9) 



P(i) = zt 



Y<k=v( w \ Coe i 'f identic^ 1 , k] D Coef ficient[C L -\ M — k — 1]\V) 



where we have defined C := E + xD. Now one can use the matrix representa- 
tion (4) to calculate (9) using MATHEMATICA. In Fig. 2 we have plotted (9) 
for two different values of q with L = 60 and M = 36. For this choice of the 
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Fig. 2. Plot of the density profile of the particles (9) on a chain of length L = 60 for 
M = 36 and two values of q (q > 1) above and below the critical point q c = 1.581. 

parameters the transition points are q c = 1.581 and q' c = 0.633. The density 
of particle has two general behaviors for q > 1. For q > q c and in the thermo- 
dynamic limit (L^oo,M^oo,p=^) the density profile of particles is a 
shock in the bulk of the chain; while in the close vicinity of the left boundary, 
it increases exponentially. The density of particles in the hight-density region 
of the shock is equal to pm g h — 1 — QT 2 - This region is separated by a rather 
sharp interface from the low-density region in which the density of particles is 
equal to pi ow = 0. The low-density region is extended over (1 — 1 _^_ 2 )L sites. 
This can be seen in Fig. 2 for q = 2.5; however, the reason that the shock 
interface is not sharp is that our calculations are not done in real thermo- 
dynamic limit. One should expect that the shock front becomes sharper and 
sharper as the length of the system L and also the number of particles on the 
chain M increase. For 1 < q < q c the density of particles in the bulk of the 
chain is constant equal to p, it drops near the right boundary exponentially 
and increases exponentially in the close vicinity of the left boundary. The ex- 
ponential behavior of the density profile of particles near the boundaries in 
this phase is due to the finiteness of the representation of the algebra (3). 
It is known that if the associated quadratic algebra of the model has finite 
dimensional representations, the density-density correlation functions cannot 
have algebraic behaviors [13]. At q — 1 one finds p(i) = p. The density profile 
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of particles in the region q < 1 is related to that of q > 1 through (7) that is 
p(i,q)=p(L-i + l,q- 1 ). (10) 



One can also study this model on a ring of length L with periodic boundary- 
conditions. Let us assume that the number of particles and therefore their 
density fluctuates and is not a constant. In this case the probability of finding 
the system in a specific configuration C = {ti, • • • ,tl} should be obtained 
from 

P(n, ■ ■ ■ , r L ) = -^trlflinD + (1 - n)E)]. (11) 



in which tr[- ■ ■] is the trace of the products of matrices. The normalization 
factor Zl can be obtained from the fact that J2c P{C) = 1- This function plays 
the role of the grand canonical partition function of the system. One can easily 
check that the quadratic algebra associated to the periodic boundary condition 
case is the same as (3) except the boundary terms which contain the vectors 
\V) and (W\; therefore, we can still use the algebra (3) and its representation 
(4) to calculate (11). As we mentioned the grand canonical partition function 
of the system can be defined as 

L 

Z l= E tr[H(T i D + (l-T i )E)]=tr[(E + D) L }. (12) 

{r=0,l} i=l 

This can easily be calculated and we find 

Z L = 1 + q~ 2L + q 2L + (1 + A) L . (13) 



The study of the zeros of (13) as a function of q in the thermodynamic limit 
L — > oo shows that they approach the positive real-q axis at two different 
points q c = a/1 + A (1 < q c < oo) and q' c = ^= A (0 < q' c < 1). Moreover, 
the zeros approach the positive real-q axis at angle |; therefore, unlike the 
reflecting boundaries case both phase transitions in this case are of first-order. 
In the steady state the density profile of the particles on the ring is flat and 
equal to p{i) = for q' c < q < q c and p{i) = for q > q c and q < q' c . One 
can take the total density of particles on the ring as an order parameter and 
since it changes discontinuously over the transition points, as the Yang-Lee 
theory predicts, the phase transitions are of first-order. One should note that 
in comparison to the reflecting boundaries case not only the geometry of the 
model is changed but also the number of particles in this case in not a con- 
served quantity. This means that we should not expect the transitions in these 
models to be similar even if we take the thermodynamic limit L — > oo. The 
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reason for the existence of two phase transition points is again the symmetry 
of the model. 

In this letter we have studied a branching-coalescing model in which parti- 
cles hop, coagulate and decoagulate on one-dimensional lattice of length L. 
We have restricted ourself to the case where the total number of particles 
on the chain is constant. For this we have worked in the canonical ensemble 
in which the number of particles M is fixed; therefore, the density of parti- 
cles on the chain p = ™ is constant. The Yang-Lee theory predicts that the 
model has two second-order phase transitions. Both phase transition points 
are determined by the density of particles on the system p. The study of the 
mean particle concentration at each site of the chain for q > 1 shows that 
the density profile of the particles has a shock-like structure in the region 
q > q c = jl_ p - The exception is near the left boundary where the density of 
particles increases exponentially. This is the first time that shocks are seen in 
one-dimensional reaction-diffusion models with reflecting boundaries. In the 
region 1 < q < q c = ^= the density profile of the particles is constant in 
the bulk of the chain; however, near the left (right) boundary it increases 
(decreases) exponentially. The exact form of the correlation lengths is under 
investigation. Our numerical investigations also show that the width of the 
shock scales as Lr v with v — |. In the thermodynamic limit L — > oo the 
shock width goes to zero and one finds a very sharp shock interface. Since 
the system is invariant under the transformation (7), the density profile of 
the particles for q < 1 can be obtained from (10). We also studied the peri- 
odic boundary case and found that the system possess two first-order phase 
transitions which are determined by the values of A for fixed value of q. In 
this case the density profile of particles is flat everywhere on the lattice and 
is either equal to or zero. The formulas (6) and (9) provide us with a 
simple and general tool for numerical study of the phase transitions and also 
the particle concentration behaviors of one-dimensional stochastic models in 
canonical ensemble for which a finite- or infinite-dimensional representation 
of the associate algebra exists. 
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